function [x1, x2, u] = td(v, r, h0, h)

len = length(v) 
x1=zeros(1,len);
x2=zeros(1,len);
u=zeros(1,len);
for k=1:len-1
    
u(k) = fhan(x1(k)-v(k), x2(k), r, h0);
x1(k+1) = x1(k) + h*x2(k);
x2(k+1) = x2(k) + h*u(k);

end